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ABSTRACT 

We present a statistical analysis of the millimeter- wavelength properties of 1.4 GHz- 
selected sources and a detection of the Sunyaev-Zel’dovich effect associated with the 
halos that host them. The Atacama Cosmology Telescope (ACT) has conducted a 
survey at 148 GHz, 218 GHz and 277 GHz along the celestial equator. Using samples 
of radio sources selected at 1.4 GHz from the Faint Images of the Radio Sky at Twenty- 
Centimeters (FIRST) Survey and the National Radio Astronomy Observatory Very 
Large Array Sky Survey (NVSS), we measure the stacked 148, 218 and 277 GHz 
flux densities for sources with 1.4 GHz flux densities ranging from 5 to 200 mJy. At 
these flux densities, the radio source population is dominated by active galactic nuclei 
(AGN), with both steep and flat spectrum populations, which have combined radio-to- 
millimeter spectral indices ranging from 0.5 to 0.95, reflecting the prevalence of steep 
spectrum sources at high flux densities and the presence of flat spectrum sources at 
lower flux densities. The thermal Sunyaev-Zel’dovich (SZ) effect associated with the 
halos that host the AGN is detected at the 5cr level through its spectral signature. 
When we compare the SZ effect with weak lensing measurements of radio galaxies, we 
find that the relation between the two is consistent with that measured by Planck for 
local bright galaxies. We present a detection of the SZ effect in some of the lowest mass 
halos (average M200 ~ 10 13 ) studied to date. This detection is particularly 

important in the context of galaxy evolution models, as it confirms that galaxies with 
radio AGN also typically support hot gaseous halos. With HerscheP observations, 
we show that the SZ detection is not significantly contaminated by dusty galaxies 
or by dust associated with the AGN or galaxies hosting the AGN. We show that 
5 mJy < Si. 4 < 200 mJy radio sources contribute £(£ + l)Cg/( 27 r) = 0.37 d= 0.03 pK 2 
to the angular power spectrum at £ = 3000 at 148 GHz, after accounting for the SZ 
effect associated with their host halos. 

Key words: galaxies: active - galaxies: statistics - galaxies: haloes - galaxies: clusters 
- radio continuum: galaxies. 


1 INTRODUCTION 

In the standard cosmological structure formation scenario, 
astrophysical structures (galaxies, galaxy groups and clus- 
ters) form when gas condenses at the centers of hierarchi- 
cally merging dark matter halos (White Sz Rees 1978). The 
accretion is characterized according to two modes: radia- 
tively efficient and radiatively inefficient (Silk 1977; Bin- 
ney 1977; Rees Sz Ostriker 1977). In the radiatively effi- 
cient rapid-cooling mode, the gas is shock-heated within the 
virial radius of the dark matter halo, and as it cools it ac- 
cretes onto the galactic disk, contributing to star formation. 
In the hydrostatic (hot) mode, the gas is shock heated be- 
fore reaching the virial radius, and a hydrostatically sup- 
ported gaseous halo forms. This gaseous halo radiates in 
X-ray wavelengths, and as it radiatively cools, it has the 
capacity to form stars. However, large cooling flows accom- 
panied by prodigious star formation are not observed: the 
predominant stellar populations at the centers of the most 
massive halos are old and red, and X-ray measurements in- 
dicate a lack of cool gas compared to the level expected 
from the cooling (e.g., Oegerle et al. 2001; Peterson et al. 
2003). One mechanism invoked to explain this problem is 
feedback from active galactic nuclei (AGN). Radiatively in- 
efficient accretion can power a “radio mode” AGN, where 
the gas that accretes onto the AGN ultimately comes from 
the hot gaseous halo. Because the radio AGN can deposit 
energy back into this gaseous halo, a feedback mechanism 
can form that regulates the gas in a way that prevents both 
catastrophic cooling and excessive star formation. Evidence 


for the existence of such feedback comes from a number of 
different observations (for a review, see McNamara Sz Nulsen 
2007). Cavities or bubbles are observed in the X-ray halos 
of nearby galaxy clusters that are filled with radio emission 
(e.g., Fabian et al. 2000; McNamara et al. 2000). Galaxy 
evolution models that include this feedback mechanism can 
successfully reproduce observed galaxy luminosity functions 
(e.g., Bower et al. 2006; Croton et al. 2006; Sijacki et al. 
2007). The correlation of the prevalence of low luminosity 
radio AGN with stellar mass, black hole mass, galaxy veloc- 
ity dispersion and color also provides support for the pres- 
ence of radio mode feedback (e.g., Best et al. 2005; Best Sz 
Heckman 2012). That radio- loud AGN tend to reside in rel- 
atively massive (~ 10 13 dark matter halos has been 

confirmed by studies of the clustering of galaxies around ra- 
dio AGN (Overzier et al. 2003; Magliocchetti et al. 2004) 
and by stacked weak lensing (Mandelbaum et al. 2009). 

Sensitive millimeter-wave observations have the poten- 
tial to reveal the Sunyaev-Zel’dovich effect of the hot gas 
associated with radio AGN. The Sunyaev-Zel’dovich (SZ) 
effect is a spectral distortion in the cosmic microwave back- 
ground (CMB) that occurs when CMB photons inverse- 
Compton scatter off of the hot ionized gas associated 
with dark matter halos. The thermal SZ effect appears 
as a temperature decrement below (and increment above) 
~220 GHz, at which frequency there is no distortion (Sun- 
yaev Sz Zel’dovich 1970). The amplitude is effectively inde- 
pendent of redshift (for a review, see Carlstrom, Holder 
Sz Reese 2002). Given the mass of the AGN-hosting dark 
matter halos, a simple extrapolation of SZ-mass relations 
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derived from more massive halos (e.g., Sifon et al. 2013) re- 
veals that the SZ effect in 1.4 GHz selected AGN should 
offset the non-thermal millimeter emission at 148 GHz by a 
significant (~50%) fraction. The Planck Collaboration et al. 
(2012) recently reported a detection of the SZ effect of halos 
of similar mass that they selected via stellar mass. These ob- 
servations suggest the feasibility of measuring the SZ effect 
from the halos that host radio-loud AGN with existing data 
sets. If the SZ effect can be separated from the synchrotron 
source spectral behavior, we can provide evidence of the ex- 
istence of hot gas atmospheres in dark matter halos hosting 
galaxies with radio loud AGN. 

At flux densities above 1 mjy, counts of 1.4 GHz ra- 
dio sources are dominated by AGN that span a wide range 
in redshift, with median redshift ~ 1 (Condon 1989). The 
synchrotron spectral behavior of radio sources is complex 
and varies significantly from source to source (e.g., de 
Zotti et al. 2010; Sadler et al. 2006; Sajina et al. 2011; Lin 
et al. 2009). Models describing radio source spectral behav- 
ior often characterize sources into two types: fiat-spectrum 
(a ~ 0, where a is defined such that S oc iz~ a ) and steep- 
spectrum (a ~ 0.8). Optically thin synchrotron emission, 
usually associated with radio lobes, results in a steep spec- 
trum source, while optically thick emission, usually asso- 
ciated with compact cores, results in a flattened spectrum 
source (e.g., de Zotti et al. 2010). Using optical observa- 
tions and radio data at multiple frequencies, the counts of 

1.4 GHz sources have been successfully modeled with radio 
luminosity functions associated with these two populations 
(e.g., Condon 1984; Peacock 1985; Danese et al. 1987; de 
Zotti et al. 2005; Massardi et al. 2011), with the flat spec- 
trum population emerging as faint sources at z = 1. These 
radio luminosity functions have been extrapolated, assum- 
ing a spectral index, to the millimeter regime where they 
have been used to predict the source population, especially 
for use in foreground models of CMB data (CMB; e.g., Tof- 
folatti et al. 1998; de Zotti et al. 2005; Sehgal et al. 2010; 
Tucci et al. 2011). In addition to the synchrotron emission 
from the AGN, the AGN host galaxies can contain dusty 
star forming regions, as have been previously investigated 
by e.g., Seymour et al. (2011). With the advent of large scale 
millimeter- wave surveys with mJy sensitivities and the avail- 
ability of IR surveys with Herschel , we can now directly test 
the spectral behavior of 1.4 GHz-selected radio loud AGN 
across two decades of frequency. 

In addition to the AGN population that dominates at 

51 . 4 > 1 mjy, deep radio surveys at 1.4 GHz (20 cm) 
(e.g., Condon 1984; Becker, White Sz Helfand 1995; Condon 
et al. 1998; Hopkins et al. 1998; Mauch et al. 2003) have 
revealed an excess number of sources at faint radio flux den- 
sities, typically with an upturn in the number counts below 
~ 1 mjy. The observed excess has been mostly attributed to 
the emergence of a population of spiral galaxies undergoing 
significant star formation (e.g., Hopkins et al. 1998; Con- 
don 1992), although with some contribution from low-power 
AGN (e.g., Seymour et al. 2008). These star- forming galaxies 
(SFGs) are the predominant population below 0.1 mjy and 
account for > 20% of the population below S' ~ 1 mjy. The 
radio emission from SFGs comes from synchrotron radia- 
tion associated with supernova remnants. In the far-infrared 
(FIR), and similarly into the submillimeter regime, SFG 
emission is predominantly thermal from the absorption and 


re-emission of starlight by dust grains and is described as a 
gray body with a spectrum rising with frequency (S ~ u 3 ' 5 ; 
Draine 2003). 

In this study, we stack 148 GHz, 218 GHz and 277 GHz 
observations from the Atacama Cosmology Telescope (ACT) 
corresponding to the locations of 1.4 GHz selected radio 
sources to determine the average millimeter emission associ- 
ated with radio sources as a function of 1.4 GHz flux density. 
We also utilize radio data from surveys conducted by Parkes 
and the Green Bank Telescope and infrared data from sur- 
veys conducted by Herschel to constrain the median spec- 
tral energy distributions (SEDs) of radio sources over the 
frequency range of 1.4 GHz up to 1200 GHz. We use Faint 
Images of the Radio Sky at Twenty-Centimeters (FIRST) 
sources matched to NRAO VLA Sky Survey (NVSS) sources 
through the cross-matched catalogs provided by Best Sz 
Heckman (2012) and Kimball Sz Ivezic (2008). This enables 
us to benefit from both the l" astrometric precision of the 
FIRST survey, and the ability of the larger (45" FWHM) 
NVSS beam to recover more extended radio emission. By 
modeling the observed flux densities, we estimate the aver- 
age spectral behavior of low-frequency sources from 1.4 GHz 
to 277 GHz and find a robust detection of the SZ effect in 
radio- loud AGN halos, particularly evident at 148 GHz. 

This paper is organized as follows. Section 2 describes 
the data from ACT, the 1.4 GHz radio source samples, the 
IR data, and additional 5 GHz radio survey data used. Sec- 
tion 3 presents the median SEDs for AGN and for SFGs from 
a relatively low redshift radio source sample constructed by 
Best Sz Heckman (2012) by matching radio sources with 
galaxies with spectroscopic measurements from SDSS, as 
well as the results of modeling different components of these 
SEDs. In Section 4 we extend this analysis to a larger sample 
of radio sources from Kimball Sz Ivezic (2008) that extends 
to higher redshift, albeit a sample that lacks optical counter- 
part identification, and split the sample into bins of 1.4 GHz 
flux density. Finally, we discuss and summarize our results 
and their implications for the SZ effect for the typical radio 
AGN halo and for the CMB power spectrum in Sections 5 
and 6. Unless specified otherwise, we assume a flat ACDM 
cosmology with Qm = 0.30, SIa = 0.70. For the Hubble con- 
stant, we adopt the notation Hq = 70/i7o kms _1 Mpc _1 , 
and E(z) describes the evolution of the Hubble parameter 
and is defined as E(z) = ^^(1 + z ) 3 + Ha. 

2 DATA 

2.1 ACT Data 

ACT is a six-meter, off-axis Gregorian telescope in the At- 
acama Desert of Chile at an altitude of 5200 meters (Swetz 
et al. 2011). The location was chosen for the atmospheric 
transparency at millimeter wavelengths and the ability to 
observe both northern and southern celestial hemispheres. 
ACT had four observing seasons between 2007 and 2010 
and surveyed both a southern celestial hemisphere region 
(S = —53.5°) and an equatorial region. ACT observed si- 
multaneously in three frequency bands centered on 148 GHz 
(2.0 mm), 218 GHz (1.4 mm) and 277 GHz (1.1 mm) with 
angular resolutions of 1.4', 1.0 7 , and 0.9 7 , respectively. Each 
band had a dedicated array of 1024 transition edge sensors. 
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This paper uses the 2009-2010 148 GHz, 2009-2010 218 GHz 
and 2010 277 GHz equatorial observations. 

Prior to all other analyses described in this paper, 
the ACT data were matched-filtered (c.f., Tegmark Sz de 
Oliveira- Costa 1998) with the ACT beam (Hasselfield et al. 
2013b, for more details see below in Section 2.1.1) in order 
to increase the signal to noise of compact sources. The fil- 
tering procedure is the same as described in Marsden et al. 
(2013). Each map for each season was calibrated separately 
(for uncertainties in overall calibration, see below in Section 
2.1.1) and filtered with the beam appropriate for that season 
and band. The resulting calibrated, filtered maps were com- 
bined into a multi-season map via a weighted average, with 
the weights set for each pixel by the number of observations. 
The ACT sensitivity varies throughout the maps according 
to the depth of coverage. The typical rms noise level is 2.2, 
3.3 and 6.5 mjy for 148, 218 and 277 GHz, respectively. 

The ACT equatorial region contains detectable levels of 
Galactic cirrus emission. In order to check whether Galac- 
tic contamination affects our results, we applied the dust 
mask used in the analysis of the angular power spectrum of 
the ACT data (Das et al. 2013). This mask was generated 
using data from the Infrared Astronomical Satellite (IRAS) 
as processed by IRIS (Miville-Deschenes Sz Lagache 2005). 
Our results are not significantly affected when we mask the 
regions with strongest dust emission. 

2.1.1 Calibration 

The uncertainty on the absolute temperature calibration of 
the 148 GHz band to WMAP is 2% (at l = 700; Sievers 
et al. 2013; Hajian et al. 2011). The relative calibration be- 
tween the 148 GHz and 218 GHz bands was done by cross- 
correlation (Das et al. 2013), and the resulting calibration 
uncertainty for the 218 GHz is 2.6%. Therefore the 218 GHz 
calibration error is correlated with the 148 GHz calibration 
error. Because the CMB survey maps for the 277 GHz ar- 
ray have not yet been calibrated to WMAP, the 277 GHz 
calibration is derived from Saturn and Uranus observations 
(Hasselfield et al. 2013b) and has a larger uncertainty of 
7%, which is uncorrelated with the absolute calibration un- 
certainties of the 148 GHz and 218 GHz data. Additionally, 
errors in the assumed instrument beam and the map-making 
can propagate to uncertainty in the recovered ACT flux den- 
sities. In propagating these additional uncertainties, we fol- 
low the procedure used by Marsden et al. (2013) and arrive 
at additional uncertainties in 148, 218, and 277 GHz flux 
densities of 1.9%, 2.6%, and 13%, respectively. The conser- 
vative 277 GHz flux density uncertainty reflects that this 
work presents the first successful reduction of these data, 
and we expect that as the reduction process for these data 
further matures, future work will have less uncertainty on 
the 277 GHz flux density scale. 

As discussed in Diinner et al. (2012), the ACT mapmak- 
ing pipeline recovers the flux densities of bright simulated 
sources to 1%. During the mapmaking procedure, bright 
sources are identified in a first-pass map and removed from 
the time-stream data before making the final map. However, 
faint sources do not have this two-step procedure applied to 
them, so in principle their flux density recovery could differ 
from that of the bright sources. We have tested the flux den- 
sity recovery for faint sources by adding simulated sources to 


the data and putting them through the mapmaking proce- 
dures, and we find that the flux densities of faint (< 10 mjy) 
sources are also recovered to 1%. 

Combining the effects of calibration, beam measure- 
ment, and map- making, the flux density uncertainties for 
the 148 GHz band and 218 GHz band are 3%, 5%, respec- 
tively, and the correlated component is 3%. The total error 
for the 277 GHz band is 15%, which is to a good approx- 
imation uncorrelated with flux density errors in the lower 
two ACT frequency bands. The full covariance is used when 
modeling flux density data in Sections 3 and 4.2. 

2.2 Radio Source Samples 

The FIRST survey mapped over 10,000 deg 2 at 1.4 GHz (20 
cm) from 1993 to 2004 using the VLA with imaging res- 
olution of 5 ;/ (Becker, White Sz Helfand 1995; White et al. 
1997). The resulting source catalog has astrometric accu- 
racy of l" . The overall limiting flux density threshold ranges 
from 0.75 mjy to 1 mjy in its equatorial region. In addi- 
tion to FIRST, the National Radio Astronomy Observatory 
(NRAO) VLA Sky Survey (NVSS; Condon 1989) was con- 
ducted at 1.4 GHz and covers the full sky above —40° dec- 
lination to a depth of ~ 2.5 mjy. 

Best Sz Heckman (2012) construct a sample of radio 
sources with spectroscopic redshifts by matching galaxies 
from the Sloan Digital Sky Survey (SDSS) with NVSS and 
FIRST. The catalog has a limiting radio flux density of 
5 mjy. The sample is split into star- forming galaxies and 
AGN according to radio and optical spectroscopic proper- 
ties: the 4000 A break strengths and the ratio of radio lu- 
minosity to stellar mass, the ratio of radio to emission line 
luminosity, and a standard ‘BPT’ emission-line diagnostic 
(Baldwin, Phillips Sz Terlevich 1981; Kauffmann et al. 2003). 
We use 667 radio-loud AGN in the region of overlap with 
the ACT survey (with median redshift of 0.30, see Figure 
1), and 149 star-forming galaxies (with median redshift of 
0.05, see Figure 1). Only 27 of the 667 AGN are identified 
by Best Sz Heckman (2012) as high excitation radio galaxies 
based on their optical emission lines, so the population can 
be modeled as dominated by low excitation radio galaxies. 
We construct and analyze the SEDs of these populations in 
Section 3. 

The Best Sz Heckman (2012) sample has the twin advan- 
tages of providing redshifts and of differentiating between 
SFGs and AGNs. It is, however, much smaller than samples 
drawn from radio surveys without optical IDs and confined 
to relatively low redshifts. To also exploit the full statistics 
of the sensitive radio surveys and to, albeit coarsely, investi- 
gate evolution with redshift, in Section 4 we study a sample 
drawn from a combined cross-matched catalog 2 of Kimball 
Sz Ivezic (2008). Because this sample is not restricted to ra- 
dio sources matched to galaxies in SDSS, it extends to higher 
redshift. From this catalog, we use all FIRST detections 
falling within the region of overlap with the ACT survey that 
are assigned at least one match to an NVSS source. Matches 
are required to lie within 30^ of one another. Should there 
be multiple FIRST or NVSS sources within the 30^ match- 
ing radius, we use only the highest ranked match as per 

2 http://www.cv.nrao.edu/ akimball/radiocat/ 
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the matching algorithm detailed in Kimball &; Ivezic (2008). 
For the stacking analysis, we use the source locations from 
FIRST, which are measured at high astrometric precision by 
the survey’s h" FWHM beam. In our subsequent modeling 
of these sources, we use their associated NVSS flux den- 
sity measurements. The larger 45 7/ FWHM beam of NVSS 
is less prone to resolving out flux density from extended 
sources and thus provides more accurate 1.4 GHz flux den- 
sities, as well as being a closer match to the ACT beam. To 
avoid stacking on lobes of extended sources that are iden- 
tified as separate sources in the FIRST catalog and more 
generally to avoid double counting sources that fall within a 
single ACT beam, we exclude sources that have neighboring 
sources within 1 ' . This leaves 9,436 sources within the over- 
lapping FIRST and ACT observing regions (~339 square 
degrees; 305° ^ a ^ 58°, —1.5° ^ 8 ^ 1.5°). We restrict the 
sample to the 4,563 sources that have 5 < Sia < 200 mjy. 
Because we model the average SED of these sources with an 
AGN model, we want to remove SFGs from the sample. We 
mask any of these sources that lie within 5' of SFGs from 
the catalog of Best Sz Heckman (2012). This eliminates 211 
sources from the sample. The expected number of SFGs for 
this sample is 140 (for comparison, 149 SFGs were identi- 
fied in Best Sz Heckman 2012), as calculated using 1.4 GHz 
count models based on luminosity functions from Dunlop Sz 
Peacock (1990) and Sadler et al. (2002). For the remaining 
sources, assumed to be radio AGN, we adopt the redshift 
distribution from de Zotti et al. (2010) (which was fit to 
data from Brookes et al. 2008) when modeling the SED. 
This distribution is shown in Figure 1 and has a median 
redshift of 1.06 (compared to 0.30 for the Best Sz Heckman 
2012 sample). 

Sources with ACT flux density greater than three times 
the rms noise level of the map in any of the three bands 
are identified and excluded from the stacking analyses. Be- 
cause the noise is determined locally, the flux density at 
which these cuts are drawn varies with position in the map, 
but typical values for the flux density threshold are 6 mjy at 
148 GHz, 9 mjy at 218 GHz, and 18 mjy at 277 GHz. These 
sources represent outlier objects that are likely bright based 
on orientation, such as blazars (Urry Sz Padovani 1995) or 
based on chance alignments, such as rare, lensed star form- 
ing galaxies (Negrello et al. 2010; Vieira et al. 2010; Marsden 
et al. 2013). By removing these millimeter bright sources, we 
reduce known inclination and lensing dependent selection ef- 
fects present in a small subset of the sources. Similarly, nega- 
tive 3cr deviations from the mean flux density in either band 
are also excluded, which will exclude large galaxy clusters 
with SZ decrements at 148 GHz and reduce the bias po- 
tentially introduced by excluding positive deviations in flux 
density from the mean. These cuts based on the ACT flux 
densities exclude 219 (4.8%) radio sources from the stack- 
ing analysis of the Kimball Sz Ivezic (2008) sample and 23 
(3.3%) AGN from the Best Sz Heckman (2012) sample. 

The synchrotron emission from radio-loud AGN is 
known to be variable. The FIRST, NVSS and ACT surveys 
were not simultaneous, so variability could affect the inferred 
spectral indices. Because the radio sources were selected 
based on their 1.4 GHz flux densities, variability would pref- 
erentially bias the 1.4 GHz flux density high relative to 
the 148 and 218 GHz flux densities, effectively steepening 
the average spectral index. However, the 1.4 GHz selected 



Figure 1. The black histogram shows the redshift distribution 
of radio AGN from Best Sz Heckman (2012), with selection as 
described in Section 2.2. The red histogram shows the redshift 
distribution of SFGs from Best Sz Heckman (2012). The black 
curve shows the redshift distribution we assume for the Kimball 
Sz Ivezic (2008) sample, which is given in de Zotti et al. (2010) 
and based on data from Brookes et al. (2008). The dashed lines 
indicate the median values for each sample. 

sources do not tend to vary as much as sources selected at 
ACT frequencies (which are more likely to be highly variable 
blazars, e.g., Marriage et al. 2011). For example, Thyagara- 
jan et al. (2011) look for variability within the FIRST data, 
using three potential criteria to determine variability: the 
distribution of the peak flux density at different times de- 
viates significantly (> 5 a) from a normal distribution, the 
maximum deviation of the peak flux density from the mean 
exceeds 5cr, and the largest variation between data points 
on a light curve exceeds 6<r. They find that for FIRST radio 
sources with counterparts identified as SDSS galaxies, the 
fraction that is variable at 1.4 GHz is 0.6%; the correspond- 
ing fraction for FIRST sources matched to SDSS quasars is 
1 %. 

2.3 Infrared data 

We investigate the ensemble SED properties of these AGN 
by calculating the median flux densities at their positions 
across a wide range of multi-wavelength data sets. To in- 
vestigate the contribution of dust to the SED of the ra- 
dio sources, we take advantage of surveys conducted by 
the Herschel Space Observatory (Pilbratt et al. 2010) using 
the Spectral and Photometric Imaging REceiver instrument 
(SPIRE; Griffin et al. 2010) that overlap with the ACT sur- 
vey region: HerMES Large-Mode Survey (HeLMS), which 
is part of the Herschel Multi-Tiered Extragalactic Survey 
(HerMES; Oliver et al. 2012), and the publicly available 
Herschel Stripe 82 Survey 3 (HerS; Viero et al. 2013a). 

SPIRE has three bands, centered at approximately 500, 
350 and 250 pm (corresponding to 600, 857 and 1200 GHz, 
respectively). The maps used are made with SANEPIC 
(Patanchon et al. 2008). Of the radio loud AGN in the ACT 
region from Best Sz Heckman (2012), 384 (58%) fall within 

3 http:/ /www. astro. caltech.edu/hers/ 
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the HerS or HerMES survey regions, as do 80 SFGs (54%). 
Of the radio sources in the sample from Kimball & Ivezic 
(2008), 2,123 fall within the HerS or HerMES survey regions. 

2.4 Flux density measurements 

At 1.4 GHz, we use the cataloged NVSS flux densities corre- 
sponding to the Best & Heckman (2012) and Kimball & 
Ivezic (2008) sources. Data from the Parkes-MIT-NRAO 
survey 4 (PMN; Griffith & Wright 1993; Tasker et al. 1994) 
and a Green Bank Telescope survey (Condon et al. 1994) 
at 4.85 GHz provide an additional constraint on the ensem- 
ble radio spectral index. For PMN and GBT as well as the 
Herschel surveys, we measure the flux densities in the map 
of each survey at the source positions, which correspond to 
the positions of the optical counterparts for Best & Heckman 
(2012) and to the positions of FIRST sources for Kimball & 
Ivezic (2008). 

For the ACT data, as discussed in Marriage et al. 
(2011), given the form of the filter, the source-centered value 
of the filtered map multiplied by the solid angle of the beam 
is the source flux density. To calculate the flux densities from 
the map, ACT pixel values are corrected by a factor that ac- 
counts for averaging of the instrument beam peak over the 
0.5 ; -square pixel of the ACT maps. An additional correction 
factor is applied to account for the fact that measured flux 
density depends on the location of a source within a pixel. 
For example, a source will have a lower measured flux den- 
sity if it is located at the junction of two pixels instead of 
the center of a pixel. We make no effort to correct this mis- 
centering effect on a source- by-source basis, as the associated 
per-source dispersion is much lower than the rms noise level 
and averages down to a negligible level in the stack. Instead 
we apply an average correction given that any source has 
equal probability of falling anywhere within the O^-square 
pixel. These factors correspond to 1.06 at 148 GHz, 1.10 at 
218 GHz, and 1.14 at 277 GHz. 

We performed null tests by calculating the weighted 
average flux densities of randomly selected, source-free lo- 
cations within the ACT+FIRST overlap region in all three 
ACT frequency bands. These were computed in the 7 flux 
density bins used in Section 4 for 1000 trials of 4,344 sam- 
ples each. These null tests were found to be consistent with 
no signal with a x 2 of 5.7, 11.3 and 3.4 at 148 GHz, 218 GHz 
and 277 GHz, respectively, each with 7 degrees of freedom. 
The corresponding probability of a random realization to 
exceed the observed 148 GHz, 218 GHz and 277 GHz x 2 
estimates are 0.57, 0.13 and 0.85 respectively. 

Details on the stacking and uncertainty estimation for 
our analysis of the Best & Heckman (2012) sample and of 
the Kimball Sz Ivezic (2008) sample are provided in Sections 
3 and 4, respectively. 


2.2), which has spectroscopic redshift measurements (me- 
dian z = 0.30) and was categorized into populations of AGN 
and SFGs. In this section, we model each of these popula- 
tions independently. 

3.1 Stacked flux density measurements 

The SED models (described in Section 3.2.1 and 3.3.1) were 
fit to the median of the flux densities (calculated as de- 
scribed in Section 2.4) for each band. The distribution of 
the flux densities of the sources (particularly at 1.4 GHz) 
is very broad, so we use the median flux densities in each 
band in order to lessen the influence of bright outliers. For 
each source in the Best & Heckman (2012) catalog, we use 
all data available. Some sources lie outside the footprints of 
the Herschel surveys but fall within the ACT survey region, 
and we include their radio and millimeter flux densities in 
this analysis. Thus the infrared median flux densities are 
drawn from a smaller sample of sources than the radio and 
millimeter median flux densities, but all sources are drawn 
from the same parent sample, with the IR sample restricted 
only by sky area, so selection differences should not be sig- 
nificant. Indeed, the radio and millimeter median flux den- 
sities of the full sample are consistent within la with the 
median flux densities of the subsample that lie within the 
Herschel survey areas. The infrared surveys were treated as 
interchangeable, as the median flux densities for radio AGN 
in HerS are within the la uncertainties of the median flux 
densities of radio AGN in HerMES for each Herschel band. 

The covariance matrix for the median flux densities was 
constructed via a bootstrap resampling of the source flux 
densities. Because the source flux densities are likely to be 
correlated between bands (a source bright in one Herschel 
band is likely to be bright in other Herschel bands), we re- 
tained off-diagonal elements of the covariance matrix. We 
have taken the overall photometric error for the SPIRE pho- 
tometer to be 7%, of which 5% is correlated between SPIRE 
bands. We added the square of this photometric uncertainty 
to the variance of the Herschel bands and the square of the 
correlated component to the covariance between Herschel 
bands, and we similarly included the uncertainties in the 
ACT calibration described in Section 2.1.1 in the appropri- 
ate diagonal and off-diagonal terms of the covariance matrix. 
This covariance matrix was used to calculate the likelihood 
of the data given the model according to the following for- 
mulation: 

— 2 lnT ex A T C 1 A (1) 

where A is equal to the data (median flux densities for each 
band) minus the model evaluated for a given set of parame- 
ters and C is the covariance matrix (which does not depend 
on the model parameters). We used an MCMC analysis to 
determine the best fit parameters and associated uncertain- 
ties for the model. 


3 SPECTRAL ENERGY DISTRIBUTION 

CONSTRUCTION AND MODELING 3 . 2 Radio i oud AGN 

In this section, we investigate the median SED for the Best & 3.2.1 AGN spectral energy distribution model 

Heckman (2012) sample (with selection described in Section 

We construct a model for the SED of the radio AGN con- 
taining components for the synchrotron, SZ effect and dust 

4 http://www.parkes.atnf.csiro.au/observing/databases/pmn/pmn.htmlontributions. The amount of free-free emission, which can 
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contribute to the millimeter regime for star forming galaxies, 
is expected to be negligible, especially given how relatively 
little flux density it contributes to the spectral energy dis- 
tributions of the much lower redshift star-forming galaxies 
in Section 3.3. The synchrotron component is modeled with 
a single spectral index: S oc . Although the spectral in- 
dex may steepen due to electron aging, no deviation in a 
is significantly detected in our analysis of the larger sample 
(Section 4.3). The frequency dependence of the SZ spec- 
tral distortion of the CMB (for units of specific intensity) 
is given by the following equation (as in Carlstrom, Holder 
Sz Reese 2002, with derivation from Sunyaev &; Zel’dovich 
1970; Sunyaev Sz Zeldovich 1972): 

s ( x ) & ( 1 + S sz(x,T e )) (2) 

where the dimensionless frequency x = . — and 

5sz(x,T e ) is the relativistic correction to the frequency de- 
pendence. The relativistic correction is typically only ap- 
preciable for massive galaxy clusters, which comprise only 
a small fraction of the host halos of the AGN. For com- 
parison, the ACT SZ-selected cluster sample analysis (Has- 
selfield et al. 2013a) uses relativistic corrections that range 
from 3 — 10% of the cluster y parameter for a sample of much 
more massive clusters (by roughly 2 orders of magnitude, as 
seen in Figure 10) than is typical for this study. 

The SZ spectral distortion described in Equation 2 is 
attributed to the hydrostatically supported ionized gaseous 
halo around the AGN. There may be an additional SZ signal 
associated with a non-radiating relativistic plasma inside of 
“cocoons” (e.g., observed as X-ray cavities in galaxy clus- 
ters) formed by the AGN during radio-mode feedback. The 
hypothesis of such a relativistic plasma is motivated by the 
fact that the minimum non-thermal pressure associated with 
the relativistic electrons sourcing the observed synchrotron 
from the cocoons is estimated to be an order of magnitude 
too small to be in hydrostatic equilibrium with the surround- 
ing X-ray emitting gas (e.g., Blanton et al. 2001; Ito et al. 
2008). The spatial extent of the cocoon ( R < 50 kpc) is sig- 
nificantly smaller than the extent of the ionized gaseous halo, 
but the cocoon’s gas pressure will likely be on average higher 
and the spectral distortion different from that of the non- 
relativistic SZ effect of the larger host halo. The formation 
of a cocoon may be modeled analytically as a point explo- 
sion (Ostriker Sz McKee 1988) using self-similar solutions to 
describe shock motion and associated gas state parameters 
(Sedov 1959). With this or other models of the cocoons, the 
spectral distortions of the CMB by the associated relativis- 
tic plasma (Wright 1979) inside the cocoons can then be 
calculated (e.g., Yamada, Sugiyama Sz Silk 1999; Platania 
et al. 2002; Pfrommer, Enfilin Sz Sarazin 2005; Chatterjee 
Sz Kosowsky 2007). These calculations, together with sim- 
ulations (e.g., Chatterjee et al. 2008; Scannapieco, Thacker 
Sz Couchman 2008; Prokhorov, Antonuccio-Delogu Sz Silk 
2010; Prokhorov et al. 2012), imply a signal in many sys- 
tems that falls significantly below that of the SZ effect from 
the bulk of the non-relativistic gas of the halo and below 
the sensitivity of the current data. Prokhorov, Antonuccio- 
Delogu Sz Silk (2010) suggest that high pressure cocoons in 
high-redshift systems could produce an SZ effect at 217 GHz 
corresponding to flux densities greater than 1 mjy. Our re- 
sults rule out the possibility that such high pressure systems 


characterize the average SED of radio loud AGN. Beyond 
this observation, we do not attempt to put constraints on 
SZ from the hypothesized non-radiating relativistic plasma 
inside cocoons and instead attribute all the SZ effect in our 
model to the non-relativistic gaseous atmospheres. As shown 
in Section 3.2.2, this model is a good fit to the data given 
their current precision. 

We parameterize our model using the measured ampli- 
tude of the SZ effect at 148 GHz, which we call Asz ■ In order 
to use the same parameter to model the 277 GHz data, we 
take into account the SZ effect frequency dependence de- 
fined according to Equation 2 as well as apply a beam (and 
therefore band) dependent correction factor that arises due 
to the effect of the filter transfer function. This correction 
factor is discussed in detail in Section 5.1. For the Best Sz 
Heckman (2012) sample, we calculate this correction factor 
to be 15% (for comparison, the effective uncertainty on the 
calibration for the 277 GHz data is 15%). 

The dust contribution is modeled by calculating the me- 
dian flux density for the dust components of the source sam- 
ple, with every source being assigned a model dust spectrum 
of the form v 13 B U (T). Each gray body spectrum corresponds 
to a single AGN, and the redshift is assigned according to 
the AGN’s optical spectroscopic redshift. Thus the flux den- 
sity observed at frequency v for a source at redshift 2 is given 
by 

q (T S Lj-R M 1 + Z)) PB vQ + Z) (' T) 

1 4*(1 + *)«&(*) f (v')P B„i (T)dv' (6> 

where the limits of the integral are 300 to 2.1 xlO 4 GHz, 
d M (z) is the comoving distance and B U (T) is the black body 
spectrum for a given dust temperature. This parameteriza- 
tion is similar to that used in Addison, Dunkley Sz Bond 
(2012), except that our model only includes a single tem- 
perature dust component. Each source is assigned a tem- 
perature of 20 K and an emissivity index of (3 = 1.8. Hard- 
castle et al. (2013) fit dust SEDs to the average Herschel 
flux densities for a similarly selected sample of radio-loud 
AGN with SDSS spectra available. They found the best-fit 
(3 = 1.8. With this (3 , for a model with a single tempera- 
ture dust component, they found the best fit T = 20.3 K, 
although they also found evidence for multiple dust temper- 
ature components contributing. We fit for a single param- 
eter describing the amplitude of every spectrum, which we 
call Lir , corresponding to the model’s infrared bolometric 
luminosity of the dust component. If there is a positive cor- 
relation between the AGN and dusty IR sources, there may 
be multiple dusty sources preferentially clustered near AGN 
that could contribute to the median flux density measured, 
potentially biasing it higher than for a single source (as dis- 
cussed for galaxies that contribute to the CIB; e.g., Dole 
et al. 2006; Marsden et al. 2009; Viero et al. 2013b). Phys- 
ically correlated structures must be near the AGN in red- 
shift space, so the spectral shape is likely to be unaffected, 
although the bolometric luminosity of the dust measured 
would become an upper limit for any single component. 

3.2.2 Results of fitting the median AGN spectral energy 
distribution 

For the radio-loud AGN, we fit for the following parameters 
of the model outlined in Section 3.2.1: the radio spectral in- 
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Figure 2. The median flux densities in NVSS, PMN/GBT, ACT 
and Herschel surveys of radio AGN from Best & Heckman (2012). 
As described in Section 3.2.1, a model with synchrotron, SZ ef- 
fect and dust components is fit to the AGN data, and the model 
parameters describing the best fit (shown here) are listed in Ta- 
ble 1. The dashed lines illustrate the synchrotron plus SZ effect 
component of the model and the dust gray body component of 
the model, and the solid line illustrates the full model evaluated 
for the best fit parameters. The error bars on the medians shown 
are the diagonal elements of the covariance matrix computed via 
bootstrap sampling (off-diagonal elements were included for the 
fitting). Some of the error bars are smaller than the size of the 
symbols. The x 2 °f the fit is 5.8, with 4 degrees of freedom. 


Table 1. Best fit parameters for median AGN SED 


a 0.55 ±0.03 

A SZ 0.45 ± 0.13 mJy 

logio^W^©) 8.7 +0.1/-0.3 

Ssync 12.2 ± 0.5 mJy 


dex (a), the amplitude of the SZ effect at 148 GHz ( Asz ), 
the typical bolometric IR luminosity of the dust component 
expressed in solar units (log 10 (Lir/Lq)), and a normaliza- 
tion parameter for the synchrotron emission. The best fit 
parameters are listed in Table 1. The uncertainties quoted 
correspond to the 68% confidence intervals returned by the 
MCMC for each parameter. Figure 2 shows the median flux 
densities for each band with the best-fit model. The % 2 for 
the best fit AGN model is 5.8, with 4 degrees of freedom. 

As evident in Figure 2, the infrared data constrain the 
dust contribution to the millimeter region of the spectra to 
lie below the detected median ACT flux densities for the 
radio AGN. There is 3 cr level evidence that the SZ effect 
contributes to the shape of the spectrum in the ACT fre- 
quencies. If we set the Asz term to 0, the % 2 is 16.0, with 

Table 2. Best fit parameters for median SFG SED 


a 0.9 +0.5/-0.4 

Ssync 5.5 +0.7/-0.6 mJy 

EM ff 943 +245/-308 cm -6 pc 

Sco 2.8 ± 0.8 mJy 

log lo (L /R /L 0 ) 9.13 +0.07/-0.09 


Figure 3. The median flux densities in NVSS, PMN/GBT, ACT 
and Herschel surveys of SFGs from Best &; Heckman (2012). De- 
scribed in Section 3.3.1, a model with synchrotron, free-free, the 
CO line for 2 : = 0.05 (rest frame frequency 230.5 GHz) and dust 
gray body emission is fit to the SFG data, and the model param- 
eters describing the best fit (shown here) are listed in Table 2. 
The dashed lines illustrate the synchrotron component, the free- 
free emission component and the dust gray body component of 
the model, and the solid line illustrates the full model evaluated 
for the best fit parameters. The error bars on the medians shown 
are the diagonal elements of the covariance matrix computed via 
bootstrap sampling (off-diagonal elements were included for the 
fitting). The error bars on the low frequency median flux densities 
are smaller than the size of the symbols. 


5 degrees of freedom. The value and significance of Asz are 
robust to changes in the assumed dust temperature (T = 
10, 15, 20, 30 K were tested) and p (P = 2.0, 1.8, 1.5, 1.0, 0.5 
were tested). 

We have also experimented with models that do not 
involve an SZ effect signal, but rather attempt to fit the 
data by altering source spectra. Unfortunately, few datasets 
constrain the spectral behavior of large samples of radio 
sources at these frequencies. The Planck mission constrains 
the millimeter-wavelength SEDs of very bright (for exam- 
ple, matched to sources with S 20 > 300 mJy) radio blazars 
(Planck Collaboration et al. 2011a) and finds some evidence 
for spectral steepening. Although we exclude such sources 
from this analysis, we nonetheless consider a model with a 
spectral index change of -0.5 at 70 GHz and allow the spec- 
tral index of the dust re-emission (p) to take on values that 
produces a minimum chi-square. For this model, p = 0.9 
produces the best fit, but the resulting x 2 is 9, compared to 
the x 2 of 6 for the preferred model (with no spectral steep- 
ening, p set to 1.8, and a term for the SZ effect amplitude). 
When p is set to 1.8, which is typical for nearby galaxies 
(Smith et al. 2013), the resulting % 2 is H- However, the 
synchrotron spectral behavior for the bright Planck sources, 
which consist mostly of blazars, differs significantly from the 
average spectral behavior of the sources in our sample. To 
allow more flexibility in the synchrotron shape, we adopt 
a model where we fix the dust spectrum p to be 1.8 and 
the SZ effect contribution to be 0, but introduce parameters 
for the location of the break in the synchrotron spectrum 
inbreak) and the amount of steepening (5a, defined such that 
a(y > inbreak) = ol + 5a). The best fit model values for a, 
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5a, and Vbreak are 0.47 4= 0.05, 0.2 +0.3/-0.5, and ~4 GHz, 
respectively. The location of the break, Vbreak, is not well 
constrained. The % 2 value is 5.8, with 3 degrees of freedom. 
If we reintroduce the SZ effect term, Asz , into the model 
and fix Vbreak to be 5 GHz, the resulting best fit parameters 
for a , 5a, and Asz are 0.484=0.04, 0.164=0.08, and 0.34=0.15, 
respectively. The x 2 is 1.3, with 3 degrees of freedom. Thus 
even with this flexibility, a model with SZ is preferred, al- 
though with much lower significance than for a model with 
a single synchrotron spectral index. We investigate this for 
the larger Kimball Sz Ivezic (2008) sample in Section 4.3, 
where we again find that the model including the SZ effect 
is preferred, and more robustly. 

3.3 SFG spectral energy distribution 

3.3.1 SFG model 

The main components of the spectra of star forming galax- 
ies include synchrotron emission that dominates at low fre- 
quencies, dust emission that dominates in the infrared, and 
free-free emission that contributes in the millimeter wave- 
length regime. We adopt the same spectral model as used 
by Peel et al. (2011), who investigate the millimeter wave- 
length SED’s of nearby late- type galaxies using data from 
Planck. The synchrotron component is modeled as S' oc iz~ a , 
and the dust component is measured as a gray body spec- 
trum with the same form as Equation 3. This model for the 
dust emission is supported by recent work: Clemens et al. 
(2013) find that single temperature (2 IK) dust spectra (me- 
dian fit /3 = 1.83) adequately describe the millimeter wave- 
length through infrared SEDs of local star- forming galaxies, 
using Planck High Frequency Instrument data in combina- 
tion with IR data from Wide-Field Infrared Survey Explorer 
( WISE ) , Spitzer, IRAS, and Herschel. We fit for the parame- 
ters a and Lir. In addition to the dust emission, the free- free 
emission is expected to contribute to the ACT bands and is 
modeled as 

Sff = 2 x 10 26 fer e e (1_ ' r s r) n^ 2 /c 2 (4) 

where Sff is expressed in Jy, k is the Boltzmann constant 
in J/K, T e is the electron temperature in K, O is the solid 
angle the galaxy subtends in steradians, tjj is the optical 
depth, v is the frequency in Hz and c is the speed of light 
in m/s. We assume T e = 8000 K, as is typical for the Milky 
Way, and calculate Q for a 10 kpc disk. The optical depth 
is calculated by 

T ff = 3.014 x 1Q~ 2 T^ 15 u~ 2 EMg gg (5) 

where v is expressed in GHz, T e in K, EMff is the emis- 
sion measure in the commonly cited units of cm - 6 pc and 
gff is the Gaunt factor, which causes curvature in the 
high frequency free-free spectrum and is approximated as 
gff = ln(4.955 x 10~ 2 /iz) 4- 1.51n(T e ) for iz in GHz and T e in 
K. The numerical factors are given by Oster (1961). We fit 
for the emission measure to characterize the contribution of 
the free-free emission to the SFG spectrum. 

In addition to these sources of continuum emission, we 
expect that the 218 GHz band will contain some contri- 
bution from the CO J(2 — 1) spectral line at 230.5 GHz. 
The redshift distribution of the SFGs spans from 0.010 to 
0.283 (see Figure 1), with median redshift of 0.047. The 


ACT 218 GHz band is 17.0 GHz wide, with central fre- 
quency at 219.7 GHz (Swetz et al. 2011). Approximating 
the band transmission as a step function, for sources in the 
range 0.01 < z < 0.09, the CO line will fall within the ACT 
band. This corresponds to 80% of the SFGs in the sample 
that fall within the ACT survey region. The true contri- 
bution of the CO line to the average SED depends on the 
detailed shape of the band transmission. In order to include 
the CO line in the SED model, we have added a parameter 
for additional flux density in the 218 GHz band. As a result, 
the 218 GHz data do not constrain the SED continuum, al- 
though obtaining a reasonable value for the CO flux density 
contributed can indicate consistency. 

3.3.2 SFG modeling results 

The SED model was fit to the median of the flux densi- 
ties of the SFGs for each band as described for the AGN 
in Section 3.2.1. We fit for the following parameters of the 
model outlined above: the radio synchrotron spectral index 
(a), a normalization parameter for the synchrotron emission 
at 1.4 GHz, the emission measure of the free- free emission 
( EMff ), the additional flux density at 218 GHz attributed 
to CO line emission ( Sco ), and the typical bolometric IR 
luminosity of the dust component expressed in solar units 
(\og 10 (LiR/ L q)). Best-fit parameters are listed in Table 2, 
and the best-fit model is shown with the median flux den- 
sities in Figure 3. The x f° r the best fit SFG model is 
4.5, with 3 degrees of freedom. As with the AGN model, 
the covariance includes both measurement uncertainty and 
intrinsic variation in the flux densities. 

We can compare the best fit values of the median SFG 
SED to published models for nearby star-forming galaxies in 
the literature. Peel et al. (2011) use Planck data to constrain 
the SEDs of M82, NGC 253 and NGC 4945. They fit for more 
parameters than our data can constrain (particularly for the 
gray body spectrum, such as the dust temperature and f3). 
Comparing with their results provides a consistency check. 
The emission measures they quote for their best fit model are 
920 4= 110 cm _6 pc for M82, 284 4= 17 cm“ 6 pc for NGC 253 
and 492 4= 81 cm _6 pc for NGC 4945. Our ensemble value 
of 943 +245/-308 cm _6 pc is in agreement with M82. Their 
best fit synchrotron spectral indices range from 1.1 4= 0.1 
to 1.6 4= 0.4, and our value of 0.9 4- 0.5/ — 0.4 also agrees 
well. Degeneracy between a and EMff is expected such that 
steeper values of a imply higher values of EMff, and this 
degeneracy is observed in the MCMC posterior distribution 
of those parameters. Including data for more bands at low 
frequency would better constrain these parameters. 

We can also compare the CO line flux density with ex- 
pectations from the literature. Bayet et al. (2006) measure 
the flux density of the CO J(2 — 1) line for nearby starburst 
NGC 253. Scaling by the square of the luminosity distances, 
the CO J( 2 — 1) flux density at the median SFG redshift 
would be 36 mjy. NGC 253 is a bright starburst, so its CO 
line flux density is likely brighter than that of typical SFGs. 
The gas density of star forming galaxies correlates with 
the star formation rate through the well-known Kennicutt- 
Schmidt relation (Schmidt 1959; Kennicutt 1998). Genzel 
et al. (2010) find a linear relation between the CO luminos- 
ity and the far IR luminosity for star forming galaxies. If we 
scale the flux density by the ratio of the bolometric lumi- 
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nosity of NGC 253 (Rice et al. 1988) to the bolometric IR 
luminosity best fit to the SFG median SED, the expected 
CO flux density implied is 3.6 mjy, which lies within la of 
the best fit CO flux density contribution to the SED. 


4 MILLIMETER WAVELENGTH BEHAVIOR 
FOR DIFFERENT RADIO SOURCE FLUX 
DENSITIES 

In this section, we investigate how the SEDs of the radio 
sources relate to their 1.4 GHz radio source flux density by 
binning the sample from Kimball & Ivezic (2008). This sam- 
ple is much larger than that used in Section 3, but it lacks 
optical counterparts and thus lacks individually measured 
spectroscopic redshifts and classification as AGN or SFGs. 
This sample also extends to higher redshift (Figure 1), en- 
abling a comparison of properties of a low redshift sample 
(discussed in Section 3) and a higher redshift sample. 

4.1 Stacked flux density measurements 

We compute the mean of the flux densities, which are calcu- 
lated as described in Section 2.4, corresponding to subsets of 
sources from the Kimball &; Ivezic (2008) radio source cat- 
alog (with selection described in Section 2.2). First, we log- 
arithmically bin sources by their associated 1.4 GHz NVSS 
flux density (Si. 4 ) into 7 bins with centers ranging from 6.4 
to 149.1 mjy. Then, for each source within a flux density 
bin, we determine the ACT flux density and the number 
of ACT observations at the FIRST position of the source, 
which determines the noise in any given pixel (see Mar- 
riage et al. 2011). Weighted by the number of observations, 
the ACT 148 GHz, 218 GHz and 277 GHz flux densities 
(S 148 , 5218 , and S 277 respectively) are then averaged to give 
the stacked ACT flux density for each 1.4 GHz flux den- 
sity bin. Weighted averages are similarly calculated for the 
GBT/PMN data and Herschel , with the weights defined as 
the inverse of the square of the uncertainties on the flux 
density measurements. For the Herschel surveys, which do 
not cover the entire region of sky used in the ACT analysis 
(2,123 out of 4,344 of the radio sources used in this analysis 
lie within Herschel survey footprints), we combine all of the 
sources across 1.4 GHz flux density bins, calculating a single 
ensemble averaged flux density for each Herschel band. Con- 
straining a single dust component also has the benefit that 
we do not need to know the redshift distribution of sources 
as a function of their 1.4 GHz flux densities and can model 
the full population as has been studied in the literature. 

The distribution of ACT stacked flux densities within 
each bin was investigated via a Monte Carlo bootstrap anal- 
ysis, resampling the flux density distribution with replace- 
ment for 1000 trials. We use the weighted averages of the ob- 
served sample flux densities, and we adopt the uncertainties 
based on the bootstrap samples. For the modeling in Sec- 
tion 4, we only include the covariance from the calibration 
uncertainties calculated in Section 2.1.1 and not bootstrap 
sampled covariance because we bin by 1.4 GHz flux den- 
sity, which limits the amplitude of the cross-frequency-band 
sample covariance relative to the dominant noise in the map 
(unlike in Section 3.2.1, which includes sources across the 
full range of 1.4 GHz flux density). 


Table 3 lists the results of stacking the PMN/GBT, 
ACT and Herschel data at the FIRST source locations 
for the radio source sample from Kimball Sz Ivezic (2008). 
Thumbnails of the stacked ACT data for each 1.4 GHz flux 
density bin of the Kimball Sz Ivezic (2008) sample are shown 
in Figure 4. Figure 5 shows the stacked ACT flux densities 
at 148 GHz, 218 GHz and 277 GHz. The stacked ACT flux 
densities associated with the FIRST sources are detected at 
^ 3 a significance at 148 and 218 GHz for all but the highest 
1.4 GHz flux density bin. The stacked ACT flux densities 
typically lie between 0.2 and 2.0 mjy. The faintest 1.4 GHz 
flux density bin contains 1,767 sources and the brightest 
contains just 61 sources. 

The average S 148 , S 218 and S 277 are always lower 
than the corresponding Si. 4 , as would be expected for 
synchrotron dominated sources with radio spectral index 
a ia -148 > 0, where flux density S oc iy~ a . However, the 
average S 218 and S 277 are greater than the average S 148 
for all but the highest corresponding bin in Si. 4 , and for 
most bins S 277 is also higher than S 2 is- In contrast, the 
synchrotron sources detected in ACT and South Pole Tele- 
scope (SPT) data, which are typically blazars and would be 
excluded from this analysis, tend to have flat spectral in- 
dices (~ 0.2) from low frequency (1.4 GHz or 5 GHz) to 
the millimeter regime, but then have falling spectral indices 
(a ~ 0.5) between the millimeter bands (Vieira et al. 2010; 
Marsden et al. 2013; Mocanu et al. 2013). There are two pos- 
sible explanations for this spectral behavior in the stacked 
ACT flux densities of radio sources: either the 1.4 GHz se- 
lected sources have, on average, an inverted millimeter spec- 
tral index a:i 48 - 2 i 8 < 0 and QU 48-277 < 0 or the SZ effect 
is causing the SED to rise from a decrement at 148 GHz 
through a null at 218 GHz to an increment at 277 GHz. We 
investigate each of these scenarios in the following sections, 
finding that the SZ effect provides a better fit to the data 
and thus is the preferred explanation. 

4.2 Spectral Energy Distribution Modeling 

We fit the model outlined in Section 3.2.1 to the weighted 
average flux densities for each 1.4 GHz flux density bin. For 
this model we allow the synchrotron spectral index a to 
vary as a function of 1.4 GHz flux density with the following 
parameterization : 

a = a 0 + /31og 10 ((5i.4>/5 p i„ ot )) (6) 

and where ao and /3 are fit as a function of Si. 4 . The central 
frequencies used here correspond to the synchrotron effec- 
tive band centers (Swetz et al. 2011). The Si . 4 normalization 
( Spivot = 149.1 mjy, corresponding to the average 1.4 GHz 
flux density for the highest bin) has been chosen to mini- 
mize degeneracy between ao and /3. Recalling that the AGN 
population can be characterized by a bimodal distribution 
including flat (a ~ 0 ) and steep (a ~ 0 . 8 ) spectrum sources 
and that the average synchrotron spectrum could also be 
affected by steepening due to electron aging, we also inves- 
tigate whether there is evidence for a break in the average 
synchrotron spectrum (see Section 4.3). We parameterize 
the Asz as constant across all 1.4 GHz flux density bins, 
motivated by previous observations that the distribution of 
radio luminosities does not correlate with host galaxy mass 
or cluster velocity dispersion (Best et al. 2005, 2007). 
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Table 3. Median flux density at radio source locations. 


Bin 

Si a Range 
(mJy) 

Nbin 

Si A 
(mJy) 

S4.8 

(mJy) 

S'l48 

(mjy) 

S218 

(mJy) 

S277 

(mJy) 

*5600 

(mJy) 

■S857 

(mJy) 

S'l200 

(mJy) 

1 

5.00 - 8.47 

1767 

6.43 

3.7 ± 0.2 

0.37 ± 0.05 

0.66 ± 0.08 

1.0 ± 0.2 




2 

8.47 - 14.3 

1092 

10.9 

5.7 ± 0.3 

0.38 ± 0.07 

0.79 ± 0.1 

1.3 ± 0.2 




3 

14.3 - 24.3 

672 

18.5 

7.8 ± 0.3 

0.59 ± 0.09 

0.83 ± 0.1 

1.5 ± 0.3 




4 

24.3 - 41.2 

412 

31.3 

12. ± 0.5 

0.64 ± 0.1 

1.1 ± 0.2 

0.96 ± 0.3 

3.9 ± 0.4 

4.4 ± 0.4 

4.4 ± 0.4 

5 

41.2 - 69.7 

222 

53.5 

21. ± 0.9 

0.94 ± 0.1 

1.4 ± 0.2 

1.5 ± 0.5 




6 

69.7 - 118. 

118 

88.4 

30. ± 1. 

1.4 ± 0.2 

1.8 ± 0.3 

2.7 ± 0.7 




7 

118. - 200. 

61 

149. 

45. ± 1. 

1.9 ± 0.3 

0.83 ± 0.4 

1.5 ± 0.9 





fThe Herschel data shown correspond to a single 1.4 GHz flux density bin containing all radio sources in the sample within 


the Herschel survey area. 
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Figure 4. Thumbnail images (0.25° x 0.25°) of filtered ACT data stacked on FIRST source locations. While the 148 GHz sources fade in 
the lower flux density bin, the 218 GHz sources remain visible. The 277 GHz map is noisier, so we have shown the 277 GHz flux densities 
stacked for the entire sample instead of splitting it into 1.4 GHz flux bins, although we evaluate the 277 GHz stacked flux densities for 
each bin for use in the analysis in Section 4.2. The 1.4 GHz flux density bins used are listed in Table 3. 


We adopt the redshift distribution based on data from 
Brookes et al. (2008) with the form given by de Zotti et al. 
(2010). The redshift distribution enters the model in the 
construction of the composite gray body spectrum (Equa- 
tion 3) and in the angular shape of the average SZ pressure 
profile, which affects how the SZ signal is sampled by the fil- 
ter as first discussed in Section 3.2.1 and further discussed in 
Section 5.1. For this assumed redshift distribution, the mea- 
sured SZ effect at 277 GHz is 12% lower than that inferred 
by the spectral behavior (Equation 2) and the SZ effect am- 
plitude at 148 GHz. This correction factor is included in the 
construction of the model, which is parameterized by the 
SZ effect amplitude at 148 GHz, Asz • Similarly, because 
the redshift distribution for this sample differs from that 
of Best Sz Heckman (2012), the same intrinsic amplitude of 
the pressure profile would produce different measurements 
for Asz due to the effect of the filtering on the pressure pro- 
files. We take this effect into account before comparing the 
integrated Y parameter in Section 5.1. 

The CO J(2 — 1) spectral line at 230.5 GHz is not likely 
to contribute to the SED because for the vast majority of 
the sources, the redshift would place the CO line outside of 
the ACT band. 

We perform an MCMC to determine the posterior 
probability distributions for the parameters oio, /3, Asz and 
l°gio {Lir/Lq), which are given uniform priors. We account 
for the full covariance of the uncertainties on the calibra- 
tion and recovery of the ACT and Herschel flux densities 
(Section 2.1.1). Figure 5 shows the model evaluated for the 
best fitting parameters along with the data. Figure 6 shows 
the resulting parameter distributions for the main param- 


eters of interest. All chains show good convergence, with 
Gelman-Rubin R - 1 parameter < 0.005. The best fit model 
has a x 2 = 31.2 with 27 degrees of freedom (PTE=0.26), so 
the model is a good fit to the data. 

The parameter estimates from this analysis shed light 
on AGN and the SZ effect from radio-loud AGN hosts. Fig- 
ure 7 shows the best fit model for a (Equation 6 ) and the 
amplitude of the SZ effect, along with the measured values 
for each 1.4 GHz flux density bin. The SZ effect associated 
with AGN hosts Asz = 0.307 it 0.052 mjy is detected with 
5 c 7 confidence. This level of SZ is consistent with expecta- 
tions given other SZ measurements of low mass systems and 
the estimated mass of the host halos of radio-active AGN. 
The data are well fit by an SZ effect term that is constant 
for all 1.4 GHz flux density bins. See Section 5 for further 
discussion of these results, including a comparison of the SZ 
effect measured for this sample with that measured for the 
Best Sz Heckman (2012) sample. 

The spectral index from 1.4 to 4.8 GHz is consistent 
with the spectral index from 1.4 to 218 GHz for sources 
across the full range of 1.4 GHz flux densities. The AGN 
parameters /3, and ao are constrained most at intermedi- 
ate to high Si. 4 , where the SZ effect represents a smaller 
fraction of the flux density (though still ~20 - 30%) than 
in fainter bins. The estimated effective AGN spectral index 
at the highest 1.4 GHz flux density bin is 0.95 ± 0.01, con- 
sistent with optically thin, “steep spectrum” synchrotron 
emission, as observed in other studies. As indicated by the 
/3 parameter, the AGN spectral index varies with 1.4 GHz 
flux density such that at the lowest 1.4 GHz flux density bin 
(Si . 4 ~ 6 mjy), a. = 0.5. This is expected for a source popu- 
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Figure 5. The PMN/GBT stacked flux densities at 4.8 GHz 
and the ACT stacked flux densities at 148 GHz, 218 GHz and 
277 GHz, with the curves corresponding to the best-fit model 
overlaid. The dashed line indicates the AGN emission and the 
solid line indicates the sum of the AGN emission and the SZ 
effect of the halos hosting the AGN ( Asz )■ All data are listed in 
Table 3. 

lation composed of some flat spectrum AGN and some steep 
spectrum AGN where the relative prevalence of one popula- 
tion over the other varies with low frequency flux density, as 
is further discussed in Section 4.3. The a = 0.55 ±0.03 mea- 
sured for the median SED of the Best Sz Heckman (2012) 
sample in Section 3.2.2 lies in the middle of the range for this 
model. However, for the distribution of NVSS flux densities 
of the Best Sz Heckman (2012) sample, the average expected 
a according to this analysis would be 0.467 ± 0.005 (statis- 
tical uncertainty from MCMC). The imposition of the re- 
quirement for optical counterparts on the Best Sz Heckman 
(2012) sample may influence the population selected such 
that the ensemble synchrotron spectra differ from the full 
radio source sample. For example, Best Sz Heckman (2012) 
preferentially exclude high redshift galaxies, whose counter- 
parts are not easily identifiable in SDSS. If there is some 
evolution in the spectra or the relative importance of flat- 
spectrum and steep-spectrum sources, then this could ac- 
count for the differences we observe between the preferred 
spectral indices in the two analyses. 

The best fitting bolometric luminosity of the radio 
sources in this sample (log 10 (Lir/Lq) = 9.96 ± 0.08) is 
significantly higher than the best fitting bolometric lumi- 



^og 10 {Li R / L 0 ) A sz [mjy] P 


Figure 6. The posterior distributions of the key model parame- 
ters, with best fit values listed. The histograms show the single- 
parameter distributions for each parameter. The dashed lines 
show the 68% confidence regions, and the solid lines show the 
95% confidence regions. The parameters correspond to the effec- 
tive AGN spectral index at the highest 1.4 GHz flux density (cko), 
how the AGN spectral index varies with 1.4 GHz flux density (/3), 
the amplitude of the SZ from ionized gas in AGN dark matter ha- 
los (Asz), an d the bolometric IR luminosity (log 10 (Lir/Lq)). 


nosity for the lower redshift Best Sz Heckman (2012) sample 
(log 10 (L/i?/Lo) = 8.7 ± 0.1/ — 0.3), as likely indicates red- 
shift evolution in the dust emission from the host galaxies. 

4.3 Flat and Steep Spectrum AGN and 
Synchrotron Spectral Steepening 

Because the sources in this work are selected at 1.4 GHz (and 
sources detected at millimeter waves are excluded), models 
describing populations of low frequency sources apply to this 
sample more directly than those describing sources at higher 
frequencies. Massardi et al. (2011) investigate source popu- 
lations from 1 — 5 GHz. As seen in previous work (de Zotti 
et al. 2010), they find that steep spectrum sources dominate 
the 1.4 GHz source counts at all flux density levels, although 
with increasing contributions from BL Lac sources (char- 
acterized by flat spectra) at low flux densities (< lOmJy) 
and from flat spectrum radio quasars at high flux densi- 
ties (> 500m Jy). Thus for this paper (which only includes 
sources below 200 mjy), we would expect to see a popula- 
tion of flat spectrum sources become more important at low 

1.4 GHz flux densities, with the caveat that the spectral be- 
havior at high frequencies may diverge from that predicted 
by the counts based on studies up to 5 GHz. This population 
would introduce a flattening of the inferred spectral index 
at low 1.4 GHz flux densities. We do indeed find that the 
population-averaged spectral index from 1.4 GHz to mil- 
limeter frequencies is flatter at low 1.4 GHz flux densities 
compared to the equivalent spectral index for sources with 
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Figure 7. Top: the effective 1.4 GHz - millimeter spectral index 
for each 1.4 GHz flux density bin (points), compared to the best- 
fit model (solid line), for which the slope and normalization are 
fit. Bottom: the amplitude of the term that models the SZ signal 
associated with AGN, compared to the best- fit model, which is 
constrained to be a constant with 1.4 GHz flux density. In both 
plots, the gray regions indicate the 1 <t uncertainty on the best- 
fit parameters of the model. For the top plot, each data point 
shown is calculated for a given 1.4 GHz flux density bin by solving 
(<5218) — (<5 i. 4)(217.6/1.4) — “ for a. The data points shown in the 
bottom plot are calculated as the difference between the measured 
and expected (5148 )• For the model fit, we use the combined 
dataset with the appropriate covariances taken into account. 


high 1.4 GHz flux densities (see Section 4.2). Conversely, 
electron aging could cause a steepening of the average syn- 
chrotron spectrum, although some models predict that only 
the brightest blazars would tend to have such steepening at 
wavelengths longer than the sub-millimeter regime (Ghis- 
ellini et al. 1998; for a list of references in support of and 
opposed to these models, see de Zotti et al. 2010). 

In order to test for the effects of spectral steepening or 
of the contribution from flat spectrum sources, we introduce 
another parameter, Ca, defined such that ot(y > Vbreak) = 
C a ot(iz < Vbreak ), where Ca can vary from —10 to 10 (and 
thus can capture a steepening, flattening, or inverted spec- 
trum) . We also introduce a parameter for the location of the 
break in the synchrotron spectrum, inbreak • The best fit pa- 
rameters for this model are shown in Figure 8 . The location 
of the break in the spectrum is not well constrained by our 
data and tends toward the ACT frequencies. Including these 
new parameters does not affect the best-fit values for any of 
the other model parameters: Asz becomes 0.302± 0.056 mjy 
(5 cr significance). When including the entire sample, from 
low 1.4 GHz flux densities to high, the best fit value for this 
new parameter is C Q = 1.06 ±0.31. The x 2 value for the fit is 



Figure 8. The posterior distributions of the key model param- 
eters for the model in which the average spectra is allowed to 
steepen, with best fit values listed. The histograms show the 
single-parameter distributions for each parameter. The dashed 
lines show the 68% confidence regions, and the solid lines show 
the 95% confidence regions. The parameters correspond to the 
amplitude of the SZ from ionized gas in AGN dark matter ha- 
los (Asz)i the amount by which the synchrotron spectrum index 
steepens ( C a ), and the frequency at which the spectrum steepens 
(i ''break)- This frequency, Vbreaki is n °t well constrained by our 
data and tends to prefer a location for the steepening above the 
lowest ACT band at 148 GHz. The amplitude of the SZ effect 
is significantly non-zero even for models in which steepening is 
allowed. 


31.0 for 25 degrees of freedom (PTE = 0.19), so the fit is not 
significantly improved by the addition of this parameter. If 
we only model the sources with Si . 4 < 10 mjy, the value for 
Ca does not change significantly. The emergence of a popu- 
lation of flat-spectrum sources could in principle contribute 
to the term in our model attributed to the SZ effect at low 
1.4 GHz flux densities, but this scenario is not preferred by 
our data when we introduce Ca. If we completely replace 
the SZ effect parameter from the model by this C a param- 
eter, the resulting x 2 = 55.3, with 26 degrees of freedom. 
Finally, we saw in Section 4.2 that for the model in which 
we include an SZ effect parameter, the average spectral in- 
dex from 1.4 GHz to 4.8 GHz agrees well with the average 
spectral index from 1.4 GHz to 218 GHz for a given 1.4 GHz 
flux density bin, implying that the simpler model without 
spectral steepening adequately describes the data. 


4.4 Galaxy Groups and Clusters 

Radio AGN occasionally reside in massive galaxy clusters 
and groups. To check whether the SZ effect inferred from the 
modeling of the data is dominated by galaxy clusters instead 
of arising from the more typical, lower-mass environments 
that host most AGN, we identify and remove radio sources 
associated with optically-selected clusters. 
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The Gaussian Mixture Brightest Cluster Galaxy (GM- 
BCG) catalog 5 (Hao et al. 2010) contains more than 55,000 
optically-selected galaxy clusters in SDSS Data Release 7 
with redshift range 0.1 < z < 0.55. While the radio source 
redshift distribution extends to higher redshift (see Figure 
1), we can estimate the contribution of clusters to the SZ ef- 
fect signal associated with radio sources based on this lower 
redshift sample. Gralla et al. (2011) found that the num- 
ber of radio sources per unit cluster mass does not evolve 
strongly with redshift. The number density of radio sources 
does evolve with redshift, implying that if anything the frac- 
tion of radio loud AGN in clusters at redshift s above z ~ 0.5 
is lower than today. Thus a significant fraction of the con- 
tribution of clusters’ SZ effect to the average SZ effect asso- 
ciated with radio AGN halos should be captured by investi- 
gating these low to intermediate redshift systems. 

Within the area of overlap between the ACT and FIRST 
surveys, there are 1,903 GMBCG clusters. Comparing this 
subset with our 1.4 GHz sample, there were 405 sources 
within a 1' projected radius of a GMBCG cluster (out of 
4,344 total sources in the sample). Of these, only 192 would 
have been included in the stacking analysis; the others would 
have been excluded because they contain multiple radio 
sources. Excluding these sources near clusters from the anal- 
ysis changes the Si 48 average values by < 5%, which is small 
compared to the errors. The model parameters (Section 4.2) 
are not significantly affected by the exclusion of sources that 
are near clusters, indicating that the Asz term is not dom- 
inated by the richest systems. 

5 DISCUSSION 

5.1 SZ effect of radio source host halos 

As presented in Sections 3.2.2 and 4.2, there is evidence 
for the SZ effect from the hot atmospheres associated with 
the halos hosting the radio-loud AGN at the 3a level for 
the sample selected from Best Sz Heckman (2012) and at 
the 5a level for the sample selected from Kimball Sz Ivezic 
(2008). In order to compare our results between samples and 
to other measurements of the SZ effect in the literature, we 
must convert the detected SZ amplitude to an integrated SZ 
signal describing the average SZ effect of the halos hosting 
the radio AGN. 

Our measured deviation in the ACT 148 GHz flux den- 
sities of radio AGN due to the SZ effect Asz is described by 
the equation 

Asz = 2-wI 0 | g{U8 GHz) | J <f>(6)y ga ,(.0)d6, (7) 

where To is the CMB intensity, g is the SZ frequency de- 
pendence (Equation 2), and <F($) is the matched filter, all 
corresponding to the 148 GHz band. Azimuthal symmetry 
of the projected SZ effect has been assumed. The function 
Ugas(0) is the dimensionless Compton parameter. The hat 
indicates it has been modified to account for the effects of 
the ACT beam and pixelization. This profile is related to the 
gas pressure profile P(r) through the line of sight projection 


5 http : / /home . f nal . gov/~ jghao/ gmbcg_sdss_catalog . html 


according to 

y g as{d) = J dsP (Vs 2 + (ife OO 0 / 05 OO) 2 ) , (8) 

where <jt is the Thomson cross section, m e c 2 is the rest 
mass energy of the electron, and R 500 and $500 refer to the 
radius and angle, respectively, within which the average den- 
sity of the dark matter halo is 500 times the critical den- 
sity, pcrit(z) = SH 2 (z)/SttG. Though we report quantities 
in this work referenced to the larger radius -R 200 , charac- 
teristic of the halo’s virial radius, we adopt the R 500 scal- 
ing for P{r) because we have chosen to model the pressure 
according to the “universal pressure profile” (UPP) of Ar- 
naud et al. (2010). The UPP is derived from X-ray obser- 
vations of clusters, for which R 500 is a more natural scale. 
Sun et al. (2011) show that an average pressure profile for 
galaxy groups agrees well with the UPP. This profile is also 
used by Planck Collaboration et al. (2012) to model the SZ 
effect of the halos hosting local bright galaxies. 

To estimate a characteristic R 500 for the radio AGN, 
we use the lensing mass estimate from Mandelbaum et al. 
(2009) (M 200 = 2.3 ±0.6 x lO 13 /^ 1 M 0 ). For R 500 and other 
relevant quantities (e.g., angular diameter distance), we cal- 
culate medians and means based on the redshift distribu- 
tions shown in Figure 1. To use the mass estimate from Man- 
delbaum et al. (2009) with the UPP Rsoo-scaled profile, we 
convert M 200 to M 500 ~ 1.5 x 10 13 hjQ Mq assuming the dark 
matter follows a Navarro, Frenk & White (1997) (NFW) 
profile, using the concentration— mass relation from Duffy 
et al. (2008). Estimated in this way, R 500 = 0.34hy O 1 Mpc 
($500 = 1.24 7 ) for the sample of Best Sz Heckman (2012), and 
R 500 = O^/i^Mpc ($500 = 0.64 ; ) for the higher redshift 
sample of Kimball Sz Ivezic (2008). The redshift distribu- 
tion of halos for which Mandelbaum et al. (2009) estimated 
a mean M 200 is best matched to the lower-redshift sample of 
Best Sz Heckman (2012). Lacking a mass estimate for radio 
AGN-hosting halos with the higher redshift distribution of 
the Kimball Sz Ivezic (2008) sample, we have simply used 
the same M 200 to estimate R 500 for both samples. 

Unfiltered brightness profiles proportional to y ga s{0) to- 
gether with the ACT 148 GHz and 277 GHz beams are 
shown in Figure 9. Since the ACT data were filtered to opti- 
mally recover point source flux, we must correct for the bias 
arising from the interaction of the extended SZ intensity pro- 
file with the matched filter. In general, the measured Asz 
will be lower than the total “SZ flux density” of the extended 
profile. To estimate the bias correction, we apply the same 
filter kernels to maps containing beam convolved, simulated 
pressure profiles, input with known amplitude, and examine 
the resulting output amplitudes. Using simulated profiles 
placed according to the redshift distributions assumed for 
each sample, we can construct the estimated distribution of 
multiplicative factor / required to correct for this bias. This 
bias factor allows us to write Asz , as in Equation 7, in terms 
of the integrated Compton Y parameter, Y 200 , defined as the 
Compton y ga s{0) integrated to $ 200 : 

Asz ~ Io | 9(148 GHz) | f^oo (9) 

We find that for the lower-redshift Best Sz Heckman (2012) 
sample, the median / corresponds to factors of 1.32 and 
1.52 for the 148 GHz and 277 GHz bands respectively. For 
the full sample of Kimball Sz Ivezic (2008) the mean corre- 
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Figure 9. The angular profiles of the projected SZ effect pres- 
sure and for the ACT beams. The dashed line shows the profile 
assumed for the sample from Best Sz Heckman (2012) (median 
z = 0.30) and the solid line shows the profile assumed for the 
sample from Kimball & Ivezic (2008) (median z = 1.06). The red 
line shows the ACT 148 GHz beam profile, and the purple line 
shows the ACT 277 GHz beam profile. 


sponding factors are 1.23 and 1.37. The ratio of the 148 GHz 
to 277 GHz bias correction factors were used in Sections 3 
(with a value of 15%) and 4 (with a value of 12%) to model 
the median and mean SZ effect, respectively. In the remain- 
der of this section, we use the bias correction at 148 GHz 
to relate the estimates of Asz for the two AGN samples to 
corresponding estimates of Y 200 • 

For the Kimball Sz Ivezic (2008) sample, these assump- 
tions yield an integrated E(z)~ 2 ^ 3 Da(z) 2 Y 2 oo = (5.44=1.2) x 
10 ~ 8 hjQ Mpc 2 where Da is the angular diameter distance. 
For the Best Sz Heckman (2012) sample, these assumptions, 
along with the redshifts measured for the sources, yield 
an integrated E(z)~ 2 ^ 3 Da{z) 2 Y 2 oo = (1.4 ±0.5) x 10~ 7 hjQ 
Mpc 2 . The integrated Y for the lower redshift Best Sz Heck- 
man (2012) sample thus exceeds that of the mostly higher 
redshift Kimball Sz Ivezic (2008) sample, although at low 
(< 2a) significance. This may indicate evolutionary growth 
in the typical halo mass and associated gaseous atmosphere 
of radio galaxies from high redshift to today. Because the the 
Compton parameter is integrated to # 200 , it follows that the 
Y 200 thus determined depends systematically on the mass 
(Y 200 oc $ 2 oo, so Y 200 oc M 200 ) due to this geometric factor. 
While the uncertainty on the mass is propagated through 
this relation to the uncertainty quoted on the measured Y 200 , 
this can also introduce a systematic uncertainty that would 
follow the same relation. Since we have assumed the same 
mass for both samples, if the high redshift halos have lower 
average mass, the underlying difference in the integrated Y 
between the high and low redshift samples would be intrin- 
sically larger than this observed difference. 

In order to characterize the level of additional system- 
atic uncertainties, we varied the assumptions made in this 


calculation about the radio source redshift distribution, the 
typical pressure profile and the typical halo concentration 
in order to test the sensitivity of the calculated integrated 

Y parameter on these assumptions. If we adopt the red- 
shift distribution of the radio sources given by Ho et al. 
(2008), the value for the Kimball Sz Ivezic (2008) sample be- 
comes E(z)~ 2 / 3 Da(z) 2 Y2oo = (4.5 ±1.0) x 10 -8 /if 0 2 Mpc 2 , 
which is in agreement with the value calculated above for the 
redshift distribution as reported by de Zotti et al. (2010). 
Adopting the halo concentration-mass relation from Neto 
et al. (2007) yields E(z)- 2/3 D a (z) 2 Y 200 = (1.5 ±0.5) x 1CT 7 

Mpc 2 for the Best & Heckman (2012) sample and 
E{z)~ 2/3 D A {z) 2 Y 2 oo = (5.7 ± 1.3) x 1CT 8 Mpc 2 for the 
Kimball Sz Ivezic (2008) sample. Finally, if we adopt the 
best-fit profile of Planck Collaboration et al. (2013), who 
measured the SZ pressure profiles of 62 massive, low redshift 
clusters, we find a resulting value of E(z)~ 2 ^ 3 Da(z) 2 Y 2 oo = 
(2.0 ±0.7) x 10 -7 hjQ Mpc 2 for the Best Sz Heckman (2012) 
sample and E{z)~ 2/3 D A (z) 2 Y 2 oo = (8.7 ± 2.0) x 10 -8 h? 2 
Mpc 2 for the Kimball Sz Ivezic (2008) sample. This pressure 
profile parameterization is more extended than the fiducial 
Arnaud et al. (2010) profile, so a given estimate of the am- 
plitude of the signal yields a larger value for the integrated 

Y parameter. 

Figure 10 shows the integrated Y parameters for both 
samples, with the low redshift Best Sz Heckman (2012) in- 
tegrated Y parameter plotted against the mass measure- 
ment for a similarly low redshift sample from a weak lensing 
analysis by Mandelbaum et al. (2009). If the true masses 
of these systems are as implied, this study has measured 
the SZ effect for some of the lowest mass halos to date, as 
most other work investigating the stacked SZ properties of 
galaxies and groups targeted ~ 10 14 M© systems (Hand et al. 
2011; Planck Collaboration et al. 2011c; Sehgal et al. 2013). 
Planck Collaboration et al. (2012) probe down to a similar 
mass regime as we do and are also consistent with extrap- 
olations from K-M scaling relations based on high (galaxy 
cluster) mass halos. While these previous studies bin ac- 
cording to mass proxy, this work computes an average value 
over what is likely a wide range of halo masses (see Section 
5.2). There is remarkable agreement among our SZ measure- 
ment for the low redshift radio galaxies, the weak lensing 
and clustering mass measurement of radio galaxies similarly 
selected, and the measurements from Planck Collaboration 
et al. (2012) for the relation between Y and mass for local 
galaxies. The higher redshift Kimball Sz Ivezic (2008) sam- 
ple integrated Y lies below the lower redshift value, possibly 
indicating evolutionary growth in the typical mass of radio 
galaxy hosts. 

5.2 Implications for galaxy formation models 

Recent observations and theory support an overall picture 
where most radio AGN are powered by radio-mode (or hy- 
drostatic mode) accretion, where hot gas from a hydrostat- 
ically supported halo ultimately cools and fuels the AGN, 
which in turn heats that gaseous halo, establishing a feed- 
back mechanism. Our measurement of the mean SZ effect 
integrated Y parameter associated with radio- AGN-hosting 
halos provides provides a direct measurement of the hot gas 
associated with this mechanism. 

Previous studies have found that the halos that host 
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Figure 10. The SZ observable versus mass. The black point cor- 
responds to the estimated integrated Y parameter that describes 
the SZ effect of halos hosting radio sources in the model that best 
fits our data for the sample from Best & Heckman (2012). The 
mass value of the black point corresponds to the masses of halos 
hosting similarly selected radio sources as determined from their 
weak lensing signal in Mandelbaum et al. (2009). The solid line 
indicates the estimated integrated Y parameter that describes 
the SZ effect of halos hosting radio sources in the model that 
best fits our data for the sample from Kimball & Ivezic (2008), 
which has a redshift distribution that extends to much higher 
redshift than the sample from Best & Heckman (2012). The gray 
region indicates the statistical uncertainty on that measurement. 
The purple data indicate the measurements for ACT SZ-detected 
galaxy clusters with dynamical mass measurements (Sifon et al. 
2013). The red data indicate results from Planck Collaboration 
et al. (2012) from measuring the SZ effect of local bright galax- 
ies through stacking analyses, and the dashed line indicates the 
scaling relation derived in that study. That the black point agrees 
well with the data at similar masses from Planck Collaboration 
et al. (2012) shows consistency between our measurements. The 
Y measurement for the higher redshift sample is lower than for 
the lower redshift sample, which may indicate evolution in the 
typical radio source host halo. 


radio-loud AGN are indeed massive, as would be implied 
by this feedback scenario. For example, Mandelbaum et al. 
(2009) measure two-point correlation functions and galaxy- 
galaxy lensing shear signals to study the dark matter halo 
masses of radio AGN selected from FIRST and NVSS and 
cross-matched with SDSS. They find that the halos hosting 
radio AGN are on average twice more massive than those 
hosting galaxies with equivalent stellar masses but without 
radio activity. Additionally, AGN without radio activity re- 
side in galaxies that, as a population, have a very different 
distribution of stellar mass (Best et al. 2005): optically se- 
lected AGN are not only present in the large stellar mass 
galaxies, but are present in galaxies with stellar mass well 
below the threshold below which we see no more radio AGN 
activity. The mean mass for halos hosting radio-loud AGN 
is (2.3 ±0.6) x 10 13 hj 0 Mq (for comparison, the mean mass 
for halos hosting optical AGN is (1.1 ±0.2) x 10 12 hjQ M©). 
Galaxy cluster BCGs are known to be preferentially radio- 


loud (Best et al. 2007; Lin &; Mohr 2007), so Mandelbaum 
et al. (2009) also repeat their analysis after excluding ra- 
dio AGN that lie within known optically selected clusters. 
While this reduces the difference between the halo masses 
of radio AGN and the control sample by ~ 15%, they still 
see evidence that the dark matter halos hosting radio AGN 
are more massive. 

Quantifying the extent to which our measurement sup- 
ports the radio mode feedback paradigm requires a predic- 
tion for the average Y parameter. In order to make a first es- 
timate based on the theoretical and observational framework 
that is currently available in the literature, we make the fol- 
lowing assumptions. We assume a stellar mass function (for 
present day halos; Cole et al. 2001), the fraction of radio- 
loud AGN as a function of stellar mass (Best et al. 2005) 
and a prescription for attributing dark matter halo mass as 
a function of stellar mass (Li et al. 2006; Mandelbaum et al. 
2009). Figure 11 illustrates the expected distribution of stel- 
lar masses (and corresponding dark matter halo masses) for 
radio AGN, optical AGN, and all galaxies, based on these 
assumptions. Although the distribution of masses for each 
population is broad, the average halo mass of radio AGN is 
significantly higher than for either optical AGN or for the 
total galaxy population. The expected average dark matter 
halo mass is M 200 ~ 2 x 10 13 M©, which is in the range 
probed by our study (Figure 10). If radio emission from jets 
were present in all active galaxies such that they would be 
selected in this sample, the average halo mass of the popu- 
lation would correspond to an undetectable (for self-similar 
scaling, Y oc M 5 ^ 3 ) SZ effect signal given the current sen- 
sitivity and sample statistics. The fact that the measured 

Y is consistent with the higher mass distribution expected 
for radio loud AGN provides a key consistency check of the 
picture that radio jets are actively providing feedback only 
to the most massive halos. Even more directly, the SZ effect 
provides evidence for the presence of ionized gas in the ha- 
los that host radio galaxies. We hope that the availability of 
this new observation will encourage the direct calculation of 

Y from simulations, thereby providing a new constraint for 
testing cosmological structure formation models with radio- 
mode AGN feedback. 

5.3 Contribution to CMB power spectrum 

Recent large millimeter surveys have enabled more sophisti- 
cated source population and spectral index modeling of high 
frequency radio sources than has previously been possible 
(e.g., Tucci et al. 2011). In a practical cosmological con- 
text, understanding the millimeter wave spectral behavior 
of radio sources is useful for modeling and removing the con- 
tribution of radio sources to high resolution measurements 
of the cosmic microwave background. Current experiments 
like Planck, SPT and ACT can use radio sources detected in 
their surveys to constrain the spectral indices of bright ra- 
dio sources (Vieira et al. 2010; Marriage et al. 2011; Planck 
Collaboration et al. 2011b; Mocanu et al. 2013). However, 
detected sources may also be easily masked from cosmo- 
logical analyses, while faint sources that fall below the de- 
tection thresholds of these experiments still contribute to 
the observed power over the angular scales they probe. A 
better understanding of the spectral behavior of such radio 
sources enables better modeling of their contribution to the 
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Figure 11. This figure illustrates the expected range of masses 
of halos hosting the radio AGN sample. Stellar and halo mass dis- 
tribution functions (per decade in mass) for all galaxies, galaxies 
with radio AGN and galaxies with optical AGN are shown by 
dotted, solid, and dashed lines, respectively. This figure is meant 
for illustration only. The stellar mass function of all galaxies (Cole 
et al. 2001) has been scaled according to the fractions for radio 
and optical AGN from Best et al. (2005) to obtain the stellar mass 
functions for the two kinds of active galaxies. The stellar mass-to- 
halo mass conversion for normal galaxies and optical AGN from 
Li et al. (2006) has been used to create the halo mass scale on the 
top x-axis. For radio AGN this should be multiplied by a factor of 
two (shifted by log 10 (2) = 0.3; Mandelbaum et al. 2009). Vertical 
lines are means of the corresponding distributions. 

power spectra and cross-power spectra for multiple frequen- 
cies (e.g., Mason et al. 2009). Characterizing the millimeter 
behavior of radio sources can also be useful for SZ surveys 
(Lin Sz Mohr 2007; Lin et al. 2009; Sayers et al. 2012; Diego 
Sz Partridge 2010), as radio sources in clusters can fill in the 
SZ decrement at frequencies below 220 GHz. 

We estimate the total power that these sources con- 
tribute to the CMB power spectrum at 148 GHz as 

n 

c[ s = ^iVi(Si - A S zf (10) 

i = 1 

where IV* is the angular number density of sources in a 
1.4 GHz flux density bin, Si is the 148 GHz flux density from 
the radio source emission for that bin, and Asz is the flux 
density associated with the SZ effect. The observed average 
flux density at 148 GHz for each bin is Si, 0 bs = Si — Asz . 
We approximate the uncertainty for each flux density bin as 
ct c ps ~ NiSf obs y/l/Ni + (2 as,i,obs/Si,obs) 2 , and the total 
uncertainty on Cf s is calculated as the sum in quadrature 
of ct c ps over all flux density bins. The values for the to- 
tal power observed (calculated using the averages for each 
flux density bin listed in Table 3) are t(t + l)Cf s /{2 -k) = 
0.14 d= 0.02 pK 2 , 0.21 ± 0.025 pK 2 and 0.6 =t 0.1 pK 2 at 
H = 3000 for 148, 218, and 277 GHz, respectively. 

As seen Equation 10, the power contributed by the 
sources at 148 GHz is partially suppressed by the SZ effect 
term, and similarly the 277 GHz power is partially enhanced 
by the SZ effect. The total power can be broken into com- 
ponents as 

Cf S = Cf xA + Cf ZxSZ + Cf xSZ (11) 


where C AxA = NiSf accounts for the power con- 

tributed by the radio sources’ emission, Cf ZxSZ ~ NtotAgz 
accounts for the power contributed by the SZ effect of the 
hosts of the radio sources, and Cf xSZ = —2 Yfi=i NiSiAsz,i 
is the cross-power. In order to calculate the Cf xA term, 
we add to the average 148 GHz flux density for each bin 
a constant corresponding to the SZ effect contribution at 
148 GHz, Asz , for the best fit model ( Asz =0.307 =t 
0.052 mjy; see Section 4.2). Thus the contribution of radio 
sources, corrected for the SZ effect of their hosts, to the CMB 
power spectrum at 148 GHz becomes t(t + l)Cf xA / (2 tt) = 
0.37 =t 0.03 pK 2 at t = 3000. This contribution is in general 
agreement with that predicted by models: for sources with 
0.3 mjy < Si 48 < 2 mjy, Tucci et al. (2011) model predic- 
tions range from 0.29 to 0.32 pK 2 and de Zotti et al. (2005) 
predict 0.35 pK 2 . Similarly, the contribution to the power 
spectrum from radio sources at 277 GHz after correcting 
for the SZ effect (including both its frequency dependence 
given in Section 3.2.1 and the effect of the filter discussed 
in Section 5.1) is l{l + l)Cf xA /(2tt) = 0.25 ± 0.07 \iK 2 
at i = 3000. The contribution from radio sources to the 
218 GHz power spectrum is unaffected by the SZ effect and 
is estimated to be t{l + l)Cf xA /(2tt) = 0.21 ± 0.025 |^K 2 
at i = 3000. For comparison, the prior on the contribution 
of radio sources to the CMB power spectrum measured by 
ACT at 148 GHz that was used by Dunkley et al. (2013) for 
sources Si 48 < 15 mjy was 2.94:0.4 pK 2 at t = 3000. The 
estimated power of sources in the SPT power spectrum after 
masking to a level of 6 mjy is 1.34=0.2 pK 2 at £ = 3000 (Re- 
ichardt et al. 2012). The amount of power that we estimate 
the radio sources in this study (5.0 mjy < Si . 4 < 200 mjy, 
0.3 mjy < Si 48 < 2 mjy) contribute after correcting for the 
SZ effect is 0.374:0.04 pK, which accounts for a fraction of 
the assumed priors (~ 25% for SPT and ~ 10% for ACT). 

We can also evaluate the power contributed by the SZ 
effect of the halos that host radio sources and the cross- 
power between the SZ effect and AGN emission compo- 
nents. Because we bin in terms of radio source power and not 
SZ effect amplitude, there is the assumption that ( A $ z ) 
{Asz) 2 • The SZ effect term becomes £{£ + l)Cf ZxSZ /(27r) = 
0.06 4: 0.004 pK 2 (statistical uncertainty) at £ = 3000. For 
comparison, Sievers et al. (2013) find that the contribution 
of the thermal SZ effect to the ACT CMB power spec- 
trum at £ = 3000 is 3.4 4= 1.4 pK 2 . The cross-power term 
is £(£ 4- l)C AxSZ /(2n) = -0.29 4: 0.07 pK 2 . 


6 CONCLUSIONS 

We have investigated the ensemble millimeter properties of 
1.4 GHz selected radio sources by stacking 148, 218, 277, 
600, 857 and 1200 GHz data from ACT and Herschel on the 
positions of radio sources selected from two joint catalogs 
of FIRST and NVSS 1.4 GHz sources. Although most of 
the radio sources are AGN, whose spectra are expected to 
fall with increasing frequency at millimeter wavelengths, we 
see evidence for an average rising spectrum. This observed 
spectral inversion in the ACT bands is attributed to the 
SZ effect of the halos that host the AGN. We construct 
SEDs for a sample of radio sources (Best Sz Heckman 2012) 
that have optical counterparts and thus known redshifts and 
classifications as either AGN or SFGs. In order to constrain 
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the contribution to the SED from dust emission, we also 
stack data from Herschel surveys, and to better constrain 
the synchrotron spectrum we include 5 GHz data from PMN 
and a GBT survey. The AGN SED is well fit by a model 
that includes synchrotron emission, the SZ effect and dust 
emission. The SFG SED is well fit by a model that includes 
synchrotron emission, free-free emission, CO line emission 
and dust emission. 

Using a larger catalog enables the construction of SEDs 
for sources binned in 1.4 GHz flux density. The radio source 
sample from Kimball & Ivezic (2008) was constructed with- 
out optical counterparts, and while spectroscopic redshifts 
are lacking, the sample includes greater numbers of sources 
(4,344 used in this study) and extends to higher redshifts. 
Even though we exclude bright ACT sources, the binned 
1.4 GHz selected radio sources are detected in the ACT 
data. The ensemble averaged millimeter synchrotron spec- 
tral index ( S oc v~ a ) for AGN is 0.95 ± 0.01 for the bright 
(~ 100 mjy) sources, with evidence for flattening toward 
fainter 1.4 GHz source flux densities. The best- fit model 
includes an SZ contribution detected at the 5 cr level. No 
relativistic SZ effect was needed to adequately describe the 
data. 

The detection of the SZ effect in the radio galaxies’ mil- 
limeter wavelength SED provides evidence that, on average, 
galaxies that host radio AGN support hot gaseous halos. 
This is important in the context of the growing support for 
the overall picture that most radio AGN are powered by 
radio-mode accretion, where AGN are fueled by the accre- 
tion of gas ultimately from the hot hydrostatically supported 
halo. In this picture, the radio jets from the AGN in turn 
provide heat that counters the cooling of the gas, shutting 
off the AGN fuel supply by establishing a feedback loop. 

We have also calculated the contribution of these ra- 
dio sources to the CMB power spectrum. At 148 GHz, 
after correcting for the SZ effect, the sources contribute 
£(£+l)C AxA /{2n) = 0.37±0.03 |J.K 2 at £ = 3000, where the 
uncertainty is dominated by the uncertainty on the contri- 
bution to the SZ effect to the measured flux densities. The 
cross-power between the SZ effect and the radio source emis- 
sion at 148 GHz is £(£ + l)C A * sz /(2tt) = -0.29i0.07 |rK 2 . 
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